Derivation of the phase field crystal model for colloidal solidification 
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The phase field crystal model is by now widely used in order to predict crystal nucleation and 
growth. For colloidal solidification with completely overdamped individual particle motion, we show 
that the phase field crystal dynamics can be derived from the microscopic Smoluchowski equation 
via dynamical density functional theory. The different underlying approximations are discussed. In 
particular, a variant of the phase field crystal model is proposed which involves less approximations 
than the standard phase field crystal model. We finally test the validity of these phase field crystal 
models against dynamical density functional theory. In particular, the velocities of a linear crystal 
front from the undercooled melt are compared as a function of the undercooling for a two-dimensional 
colloidal suspension of parallel dipoles. Good agreement is only obtained by a drastic scaling of the 
free energies in the phase field crystal model in order to match the bulk freezing transition point. 

PACS numbers: 82.70.Dd, 64.70.D-, 81.10.-h, 81.10.Aj 



I. INTRODUCTION 

Crystal growth processes are relevant for a variety of 
different problems ranging from crystallization of pro- 
teins ill] and other biological macromolecules [H over 
the construction of photonic crystals with an optical 
bandgap Q to apphcations in, e.g., metallurgy A 
full microscopic understanding of crystal growth with the 
interparticle interactions and the thermodynamic bound- 
ary conditions as the only input is still a great challenge 
since it requires a microscopic theory of freezing. Sig- 
nificant progress has been made by using the so called 
phase field crystal (PFC) method in which the tradi- 
tional phase field theory Q is generalized to a situa- 
tion with a crystal order parameter. The PFC model 
was first developed by Elder and coworkers Q and then 
subsequently applied to many other situations like inter- 
faces polycrystalline pattern formation [1, Q, crys- 
tal nucleation |10l] , commensurate- incommensurate tran- 
sitions and edge dislocations [l^ . 

On the other hand, classical density functional theory 
(DPT), which provides a microscopic theory for freezing 
in equilibrium [H, [li, [H [M [13, E d , was general- 
ized to nonequilibrium situations for colloidal particles 
with Brownian dynamics. The so called dynamical den- 
sity functional theory (DDFT) can be derived from the 
basic Smoluchowski equation and was found to be 
in very good agreement with Brownian dynamics com- 
puter simulations for dynamics in inhomogeneous situa- 
tions [20', mi, [Hi including crystal growth ^21, [11]. The 
PFC modeling possesses a static part invoking a free en- 
ergy expression for order parameters that describe the 
equilibrium bulk crystallization transition and a dynami- 
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cal part that describes the evolution of the order param- 
eters by a generalized diffusion or continuity equation. 

Recently, Elder and coworkers have derived the static 
free energy input to the PFC model from microscopic 
equilibrium densit y fu nctional theory using a truncated 
density expansion [2^ . This was supplemented by a gra- 
dient expansion in the order parameter (26l.[27j. However, 
a microscopic justification and derivation of the dynam- 
ical part is still missing. In this paper, we close this 
gap for colloidal suspensions where the individual dy- 
namics is overdamped Brownian motion. Here, the ap- 
propriate microscopic starting point is the Smoluchowski 
equation [H, [2§] from which by an adiabatic approxi- 
mation the dynamical density functional theory can be 
derived [20*]. In this paper we show that the dynam- 
ical density functional theory provides a firm theoreti- 
cal basis to derive the phase field crystal model and to 
discuss the approximations involved. We end up with 
two versions of the phase field crystal model, referred to 
as the PFCl and PFC2 models, which can both be im- 
plemented numerically with the same effort. We argue 
that the PFCl model involves less approximations than 
the standard PFC2 model, which is traditionally used in 
phase field crystal calculations. 

Finally we compare the dynamics of freezing for dy- 
namical density functional theory (DDFT) and phase 
field crystal (PFC) modeling. The system considered 
here are two-dimensional dipoles with parallel dipole mo- 
ments pointing out of their confining plane. This sys- 
tem is realized for superparamagnetic colloidal particles 
which are confined to a two-dimensional air-water inter- 
face at a pending water droplet and exposed to an exter- 
nal magnetic field [10, [3l| . This system is characterized 
by the pairwise interaction potential u{r) ~ u^/r^ , where 
uo is a parameter with dimensions of energy x volume. 
For the specific realization of two-dimensional paramag- 
netic colloids of susceptibility x exposed to a perpendicu- 
lar magnetic field B, we have uq = (xB)^/2 in Gaussian 
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units. As for all power-law interactions, the thermody- 
namics and structure depend only on one dimensionless 
coupling parameter T = uqp'^^^ /ksT, where p is the av- 
erage one-particle density and fcsT is the thermal energy. 

For our comparison between the DDFT and the two 
PFC models we consider crystal growth out of an un- 
dercooled melt into a two-dimensional triangular crystal. 
The growth velocities of a linear crystal front, which is 
cut out of a perfect hexagonal lattice, are calculated as a 
function of the undercooling for both the PFC model and 
the DDFT. After renormalization of the excess part of 
the free energy good agreement is obtained between both 
approaches implying that the PFC approach — with ap- 
propriately chosen input parameters — provides a reason- 
able and justified description framework of crystal growth 
phenomena. 

The outline of the paper is as follows: In Section [Til 
we derive the DDFT from the Smoluchowski equation 
of overdamped Brownian motion, in the same fashion as 
presented by Archer and Evans |20|. The approximate 
free energy functional by Ramakrishnan and Yussouff [33| 
is incorporated into the DDFT in Section IIIII Subse- 
quently, in Section lTVl the two versions of the PFC model 
(PFCl and PFC2) are derived from the approximate 
form of the DDFT. The different theories (the DDFT, the 
PFCl, and PFC2 models) are applied to crystal growth 
of dipoles in two dimensions in Section IVIl including the 
presentation of the system's equilibrium phase diagram 
(Subsection IVI A|) and of the non-equilibrium problem 
(the setup description in Subsection IVIBi and the re- 
sults in Subsection ETCJ. In Section IVlTl we summarize 
and conclude. 



II. DYNAMICAL DENSITY FUNCTIONAL 
THEORY (DDFT) 

We consider the overdamped dynamics of a set of N 
identical, spherical, colloidal particles, immersed in a sol- 
vent, which serves for damping and as a heat bath. As- 
suming that the particles do not interact via hydrody- 
namic forces, the N coupled Langevin equations of mo- 
tion [28j are given by 



r, = (F, + f,) 



l,...,iV, 



(1) 



where the dot denotes a time derivative and 7 = Snrjocr is 
the friction coefficient for a colloidal sphere with diameter 
cr in a fluid of viscosity 779. For particles in an external 
field V{ri, t), which interact with each other via pairwise 
additive potentials u{\ri — rj\), the deterministic force 
acting on particle i is given by 



F,({r},t) 



1 u{\r,~r,\) + V{r,,t) 



(2) 



where we denote the positions of all particles by {r} = 
{ri, . . . , rjv}. The Gaussian white noise random forces 



originating from the solvent are characterized by the first 
two moments of their distribution function, 

(Mt)) = (3) 

{lUt)fjp{t')) = 2-fkBTS^jSa.pd{t - t') , (4) 

where kgT is the thermal energy. The angle brackets de- 
note a noise average, and Greek indices indicate a com- 
ponent of the cartesian vector. Eqs. ([3]) and ([4]) fulfill 
the well-known Einstein fluctuation-dissipation relation 
yielding a short-time diffusion constant D = ksT/^. 
The set of coupled, stochastic differential equations (H]) 
for the particle coordinates corresponds to a determinis- 
tic Fokkcr-Planck equation for the A^-particle probability 
density W{{r},t) [Hill, 

Wi{r},t)^CsW{{r},t), (5) 
>Cs=7"'I]V, •[fcBrV,-F,({r},t)], (6) 

i 

which determines the probability to find the set of N par- 
ticles within a small volume around the positions {r} at 
time t, given a normalized, initial distribution W^({r}, t = 
0). The sum runs over all particles i = 1, . . . , iV. The 
continuity equation ([5]) is referred to as Smoluchowski 
equation [28j . 

For dense, strongly interacting fluids, one is typically 
not interested in the position of all individual particles 
but rather in the probability to find any particle at a 
certain vector r at time t. We therefore introduce the 
time-dependent one- and two-particle densities 

pir,t)=J2{Sir-r,m, (7) 

i 

p(2'(r,r',i)= E {S{v-v,{t))Siv-r,m, (8) 

where we dropped the superscript "(1)" on the one- 
particle density. Generally, the n-particle density is equal 
to the {N — n)-times integrated probability density W, 



,t) 



Nl 



{N -ny. 



vWi{r},t). (9) 



A deterministic equation of motion for the time evolu- 
tion of p{r, t) — Eq. pO|) — can on the one hand be di- 
rectly derived from the Langevin equations, Eq. ([1]), 
via a coordinate transformation r.; — s- /5(r, i), where 
p{r,t) = {r{t) — ri{t)) is the one-particle density 

operator, and a subsequent noise-average. This way 
was followed by Marini Bettolo Marconi and Tarazona 
(MT) based on an earlier approach by Dean 3j|. 
On the other hand, the same equation is obtained by 
integrating the Smoluchowski equation ([5]) over the po- 
sitions of iV — 1 of the TV particles and making use of 
Eqs. (O, ([5]), and ©. The latter approach was adopted 
by Archer and Evans ^2Q] . The continuity equation for 
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p{r,t) reads 

p(r,t) = • [fci3rVp(r,t) +p(r,t)VF(r,i) 

+ / dr'p^^\r,r',t)Vu{\r-r'\) 



(10) 



For noninteracting particles in zero external field, this 
equation reduces to Pick's diffusion equation. Also with 
an external field applied, Eq. (fTO|) is exactly solvable. In 
the interesting case of interacting particles, however, an 
expression for the time-dependent two-particle density 
p^^^ (r, r', t) is still needed. 

Within DDFT, p(^^(r, r',i) is approximated by a yet 

(2) 

unspecified equilibrium two-particle density pj, (r,r'); 
the latter is evaluated at a corresponding equilibrium 
fluid, in which the equilibrium density po{r) is equal 
to the instantaneous one-particle density p{r, t) of the 
nonequilibrium system. The approximation of replac- 
ing a time-dependent, nonequilibrium by an equilibrium 
correlation function, is referred to as adiabatic approx- 
imation; it goes back to Enskog [s^l, who applied it to 
the time evolution of the single-particle distribution func- 
tion in a dense gas of hard spheres. In order to render 
the instantaneous density p{r, t) an equilibrium density, 
an appropriate external potential v(y) must be applied. 
That such a potential exists for any physical density field 
p{v, t), and that it is, further, a unique functional of the 
density po(r), is stated and proved in one of the basic 
theorems of density functional theory [ssj . 

The connection to classical density functional theory 
is now made by identifying the three different terms in 
the bracket on the right-hand side of Eq. pH]) with terms 
of the form p{r)V S Fi[p] / 6 p{r) , where the Fi[p] arc differ- 
ent contributions to the Helmholtz free energy functional 
F[p{r)], which is provided by classical density functional 
theory. The functional F[p(r)] is a unique functional of 
the static one-particle density p(r) [331 . If F[p(r)] is 
known exactly it is minimized by the equilibrium one- 
particle density p(r) = Poi^), where it takes the value of 
the Helmholtz free energy F = F[po{r)]. The functional 
is divided into three terms: 

F [p{r)] = Fid [p(r)] + F,^ [p{r)] + F,^t [p{r)] . (11) 

The ideal gas part, which is of completely entropic nature 
and which yields the (first) diffusion term in Eq. pU|) . is 

Fid [p(r)] = ksT J drp(r) {in [p(r)A'^] - l} , (12) 

with A denoting the thermal de Broglie wavelength and d 
the spatial dimension. The external part corresponding 
to the second term in Eq. pU]) is given by 



Fe,t [p(r)] = J drp(r)y(r,i). 



(13) 



Finally, the excess part Fcx[p(r)], originating from the 
correlations between the particles, is generally unknown 



and must be approximated; a specific approximation is 
introduced further down. Note that the excess part is not 
the potential energy of interaction but a contribution to 
the free energy. The connection of the excess part and 
the third term on the right-hand side of Eq. pO|) is made 
through the sum rule 

-po(r)VcW(r) = {kBT)-' J dr'p^') (r, r')Vu (|r - r'|) , 

(14) 

("21 

which connects the two-particle density Pq (r,r') with 

the effective one-body potential kBTc^^\v). The lat- 
ter, in turn, is — up to a minus sign — equal to the first 
functional derivative of the excess free energy functional 
-Pox[po(r)] with respect to density. 



kBlCo (rj- ^^^^^ 



(15) 



Using Eqs. and we can therefore rewrite 

Eq. Ho]) as 



P(r,<) 



kBTW^p{r,t) + W ■ [p{r,t)VV{r,t)] 



-V • 



5p{v,t) 



(16) 



which, making use of Eqs. (fTT|) . (fT2| . and (fT3|) . reads 
a compact form 



p(r,i)V 



5F [p(r,t)] 
<5p(r,t) 



(17) 



Eq. (|17p constitutes the fundamental, nonlinear, de- 
terministic equation for the time-evolution of the one- 
particle density p(r, t) and will be referred to as DDFT 
equation henceforth. For time-independent external po- 
tentials V{y), the DDFT describes the relaxation dynam- 
ics of the density field towards equilibrium at the mini- 
mum of the Helmholtz free energy functional F[po], given 
an exact canonical excess free energy functional Fc^[p]. 
The path in the space of density fields is in general not 
the one of steepest descent, but is governed by the mass 
conservation constraint in Eq. (flT)) plj . 

Eq. p7p is a deterministic equation; it has no extra 
noise-term, as all the fluctuations — given that fex[p] is 
exact — are already taken into account. As discussed at 
large by MT [21I, |3^, the addition of a noise term in 
Eq. leads to an overcounting of fluctuations. Archer 
and Rauscher have discussed the possibility of including 
a noise term if p(r, t) is not interpreted as the ensem- 
ble averaged but as a coarse-grained (in time) probabil- 
ity distribution 'p{r,t) [s^- However, this also requires 
a replacement of the functional F[p] by a functional of 
the coarse-grained probability density, which is not the 
Helmholtz free energy functional of density functional 
theory. We will come back to this point in Section fVl 

The DDFT equation had been su gges ted earlier on 
phenomenological grounds by Evans j35| and later by 
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Dieterich et al. [3^- The same dynamical equation with 
the excess free energy functional of Ramakrishnan and 
Yussouff fsp. 39I (cf. Section [TTT| had been derived by 
Munakata " l4lj and extended to non-spherical parti- 
cles in the context of solvation dynamics by Calef and 
Wolynes [13], which was later reformulated by Chan- 
dra and Bagchi these equations are referred 
to as Smoluchowski-Vlasov or nonlinear diffusion equa- 
tions jllj , as they are derived from a Vlasov equation ]45| 
with a Fokker-Planck collision operator [4l[. However, 
MT were the first to derive the theory from the micro- 
scopic equations of motion and to make clear the con- 
tact to static DFT [21], [sl] . Similar attempts had been 
made by other authors before (K irkp atrick et al. [46l |. 
Dean j33|, and Kawasaki et al. [43,|48|), which, however, 
do not distinguish the average density p with the density 
operator p, which in turn leads to an additional noise 
term on the right-hand side of Eq. (|17p and therefore to 
an overcounting of fluctuations, given an accurate func- 
tional F[p] (see also the discussions by MT [2l|, Archer 
and Rauscher [s^l, and Lowen [iot). 

The DDFT is an approximate theory in several re- 
spects: the first and most fundamental approximation 
is the already introduced assumption of adiabatic relax- 
ation dynamics. In practice, this approximation is most 
severe in dynamical processes that are fast compared to 
the diffusive time scale of the system. To our knowl- 
edge, this issue has been studied systematically to date 
only for weak perturbations of a hard-rod fluid in one 
dimension by Penna and Tarazona [soj . The use of ap- 
proximate free energy functionals is the second funda- 
mental approximation turning out to be severe in many 
applications (cf. the next section). Third, we did only 
consider systems in which hydrodynamic interactions be- 
tween the particles play no role. The latter assumption 
can be approximately tackled by allowing for density- 
dependent friction constants 7 [5l|, which is appropri- 
ate for long-wavelength fluctuations of the density field, 
or by taking hydrodynamic interactions on the Rotne- 
Prager (two-particle) level into account, as was recently 
demonstrated by Rex and Lowen [s^ . 



III. APPROXIMATE DENSITY FUNCTIONALS 
IN THE DDFT 



For most problems including those invol ving freezing, 
Fcx[p(r)] is only known approximately [T^. l35l|: in prac- 
tice, this restriction most often constitutes the more se- 
vere approximation as compared to the adiabatic ap- 
proximation. In this paper, we follow the approach of 
Ramakrishnan and Yussouff (RY) [s^l as laid out for 
the dipolar system in 2D in reference [s^ and as al- 
ready exploited for the DDFT of the same model in ref- 
erences [2^ nil- Within the RY-approach, Fox [p(r)] is 
expanded up to second order in terms of density differ- 
ence Ap = p(r) — p around a reference fluid, where the 
fluid density p is chosen the average density of the inho- 



mogeneous system: 



' ' drdr' Ap(r) Ap(r')4'^ (r - r'; p) . (18) 



(2.) 

Here Fex(p) and Cg "^(r; p) are the excess free energy and 
the direct correlation function of the reference fluid of 
density p, respectively Despite the ease of imple- 

mentation the RY-excess free energy functional is used 
here because it will lead directly to the PFC model in the 
next section. Within the RY-approximation the DDFT 
equation now reads 



V2p(r, t) + {kBT)-^V ■ [p(r, t)VV{Y, t)] 



-V 



p(r,t)V / drV(r' 



(19) 



with D = ksT/^ the diffusion constant. 

Apart from leading to quantitatively wrong equilib- 
rium density fields and free energies the approximate den- 
sity functional might display more than one local min- 
imum, in which the system might get trapped. E.g., 
in conjunction with freezing a fiat and constant density 
profile, corresponding to the fluid state, is "metastable" 
at all temperatures for any known approximate density 
functional; starting from the fluid state, within DDFT, 
the system therefore never reaches the stable crystalline 
state which is represented through a periodically modu- 
lated density field. As discussed at length by MT [Uii], 
this failure made some groups add a noise term to 
Eq. p7|) . which is at best justified a posteriori. In par- 
ticular, it leads to the already mentioned overcounting of 
fiuctuations (see also Section IV)) . 



IV. THE PHASE FIELD CRYSTAL (PFC) 
MODEL 

As the DDFT, the PFC model is based on a free en- 
ergy functional J-"['i/'(r, t)] of a phase field '0(r, i) and 
a dynamical equation for the phase field's time evolu- 
tion similar to the DDFT equation. The PFC model 
was introduced as a phenomenological theory by Elder et 
al. d,!!!]. If the yet to be specified functional J^[^{y, t)] is 
set equal to a particular approximation of the Helmholtz 
free energy functional from density functional theory, i.e., 
J-[ip{Y, t)] = F\p{r, t) = ^p{r, t)], as was also suggested by 
Elder et al. [25] , the phase field is consequently to be in- 
terpreted as the density field, i.e., ip{r,t) = p{r,t). We 
show in this section that the commonly used PFC equa- 
tion of motion [Eq. (^7)) ] can be regarded as a par- 
ticularly simplified and further approximated version of 
the DDFT with the RY approximation to the excess free 
energy functional [cf. Eq. (fT9|) ]: consequently, the PFC 
model is derived here from the basic Langevin equations 
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of motion ([T]) for the case of overdamped dynamics. Of 
course, this reasoning only holds if the phase field ?/'(r, t) 
in the PFC model is regarded as the density field p(r, t) 
of Eq. (O, which will be assumed in this section and is 
believed to be assumed in many other papers whenever 
is set equal to F[p]. Different interpretations of 
i/j{r, t) are discussed in the next section. 

Apart from the adiabatic and the RY approximation 
the derivation goes via three further approximations; 
First, the RY excess free energy functional, Eq. (fT8|) . is 
approximated by a (local) gradient expansion. Second, 
the mobility in the dynamical equation (jl9[) . is set to 
be constant, i.e., j~^p(r,t) « 7~^p with p the average 
density of the system. Third, the ideal gas part of the 
free energy, Eq. (I12p . is approximated by its truncated 
Taylor series. According to the (additional) approxima- 
tions of the PFC model two different equations of motion 
are put forward referred to as the PFCl model, which is 
obtained after the first approximation, and as the PFC2 
model, which is obtained after the second and third ap- 
proximations. Obviously, the PFC2 model constitutes an 
approximate form of the PFCl model. 

Following the procedure suggested by Elder and 

coworkers [25], the RY excess free energy functional, 
Eq. (US]), is approximated by its gradient expansion 
within both PFC approaches: 



drAp(r) Co-CsV^ 



Ap(r) . 



(20) 



Consequently, J^cx is local in the density field, which ren- 
ders the yet to be introduced PFC equations of motion 
[approximate forms of the DDFT equation p9)l ] com- 
putationally faster to solve. The gradient expansion is 
equivalent to a Taylor-expansion of the Fourier transform 
Cq (k; p) of the two-particle direct correlation function 
introduced in Eq. (llSp . 



Co + C2k^ + C^k^ + 



(21) 



Due to rotational symmetry of the pair correlation func- 
tion the expansion is only in even powers of k. Truncating 
the expansion at fourth order, the time evolution given 
by Eq. p?]) now reads 



p(r, t) = DV^p{r, t)+DV-{ p(r, t)V (fc^T)- V(r, t) 



-(Co-C2^^+C4^^)pir,t) 



(22) 



This equation, which we refer to as PFCl model, approx- 
imates the integro-differential equation of the DDFT, 
Eq. pT|) , by a local partial differential equation of sixth 
order. Further down, we will advocate the use of this 



equation rather than of the more approximate equation 
of the PFC2 model. 

The (second) constant mobility approximation, an ad 
hoc assumption of a constant density p{r, t) = p in front 
of the functional derivative in Eq. pT]) , leads to the equa- 
tion 



ST[p{r,t)] 



(23) 



5pir,t) 

where the total free energy functional is given by 

T[p] = Fid [p] + Fext [p] + ^ox [p] ■ (24) 

Concurrently, insertion of the ideal gas part of the 
Helmholtz free energy functional, Eq. leads to a 

term, which is logarithmic in the density field. 



SF,i [p{r,t)] 
5p{v,t) 



= kBTVHn [p(r,t)A'^] 



(25) 



This term replaces the simpler diffusion term, V^p(r,t), 
in the DDFT equation (and in the PFCl model). Within 
the PFC2 model the logarithm is expanded in a power 
series about the constant density p, i.e.. 



Fid [p(r)] « ksTp J dr{i0(r,<)2 - ^-^(r,*)' 

} 



(26) 



const 



with 0(r,t) = [p{r,t) — p]/p the dimensionless density 
modulation. This leads to the standard form of the PFC 
model used in the literature. 



(r,t)- i<^(r,i)2 + i0(r,t)3 



^(fcsT)- V(r, t)-p (Co - + C4V 



,t) 



(27) 



which is henceforth referred to as the constitutive equa- 
tion of the PFC2 model. Note that the second and 
third term on the right-hand side only appear due to 
the constant-mobility assumption and are not present in 
the less approximate Eq. ((22|) . 

We will use the (standard) PFC2 model, Eq. jST]), as 
well as the more accurate PFCl model, Eq. (j^^ . and 
compare them with the DDFT, employing the RY ap- 
proximation, Eq. (|19p . Therefore we need to parametrize 
Co, 6*2, and C4 in the PFCl and PFC2 models according 
to c(k; r) in Eq. (j2ip . A particular parametrization, mo- 
tivated from the one-mode approximation to the PFC2 
model [11], is chosen and presented in SectionlVTl Before 
we come to the comparison we comment shortly on the 
use of an additional noise term in Eqs. (|^^ and (ITT)) in 
the following section. 
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V. NOISE TERM IN THE PFC EQUATION 

Typically, the PFC2 equation ((27)) is supplemented by 
a non-multiplicative Gaussian noise term ?7(r, t) 6] ful- 
filling 

(77(r,t)) = 0, (28) 
(?7(r,t)?7(r',i')> = 2kBTj-^S{r ~ r')S{t - t') . {29) 

As already pointed out in Section |TT1 such a noise term 
can not be derived in the context of DDFT, since the 
density field p(r, t) is an ensemble averaged quantity. In- 
stead, the addition of a noise term in Eq. (|16p leads to 
an overcounting of fluctuations [2l| . Therefore we point 
out that — at least for colloidal dynamics — the addition 
of a noise term is not well justified. 

However, Archer and Rauscher (STj argue on a phe- 
nomenological basis that a noise term fulfilling Eqs. pS)) 
and (|29p can be introduced if the phase field is not under- 
stood as the ensemble averaged probability density p(r, t) 
but as a coarse-grained (in time) density 

p(r,i) = / dt'K{t-t')piv,t'), (30) 

J — QO 

with K{t) a coarse-graining function of width t = 
J^dttK{t) and p{r,t) = J2^HHt) -^t)) the den- 
sity operator, as already introduced above. However, if 
this approach is followed, the functional F[p] needs to 
be replaced by a functional of the coarse-grained prob- 
ability density p(r,t), which is not the Helmholtz free 
energy functional of density functional theory and which 
is generally unknown. Second, the temperature enter- 
ing Eq. (j29p must be renormalized by a factor \/to/t 
accounting for the ratio of the microscopic time-scale tq 
and the coarse-graining time scale t. In fact, the depen- 
dence of Eq. (|29p on tq points to a conflict with a funda- 
mental assumptions of Brownian dynamics, namely that 
To is much smaller than any relevant time scale. If this 
reasoning is followed nevertheless, the non-multiplicative 
nature of the noise-term, Eq. ((28)) . comes about only after 
an approximation similar to the one of constant mobility 
in the PFC2 model whereas the less approximate PFCl 
model should be appended by a multiplicative noise term 
of the form 



77(r,i)--V^^p(r,i)77(r,t), (31) 

where the gradient assures mass conservation. Due to the 
gradient the Ito and Stratonovich calculus [11] are equiv- 
alent. For a phenomenological motivation of the accord- 
ing dynamical equation, we refer the reader to Archer 
and Rauscher |37|. In the following, we do not consider 
a fluctuating density field but restrict our study to the 
application of the deterministic DDFT, PFCl, and PFC2 
equations, respectively. Finally, we remark that a noise 
term in the case of molecular dynamics was recently dis- 
cussed by Tupper and Grant [56l |. 
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FIG. 1: (Color online) The Fourier transform c{k) of the two- 
particle direct correlation function for different coupling con- 
stants r = 10,30,50, plotted against k/p^^^. 



VI. APPLICATION OF THE DDFT AND THE 
PFC MODEL TO PROPAGATING CRYSTAL 
FRONTS 

In this and the following section, we compare the 
three different approaches to the non-equilibrium, Brow- 
nian dynamics discussed, the DDFT and the two PFC 
models, to the problem of propagating crystal fronts. 
For a better understanding of what drives the crystal 
growth, we present, first, the equilibrium phase diagram 
obtained from the underlying free energy functionals in 
Subsection I VI A) before studying the dynamics in Subsec- 
tion [VlBj 



A. The equilibrium state 

Input to the static free energy functionals, Eq. ((TT)) and 
Eq. ()24p . and, concurrently, to the dynamical theories is 
the direct correlation function of the fluid c^\r) \Et\ . 
which has been obtained for a large range of coupling 
constants < F < 62.5 from liquid-state i nteg ral equa- 
tion theory as described in references [H, Issf . In par- 

(2) 

ticular, Cq (r) has been obtained by iteratively solving 
the coupled Ornstcin-Zernicke equation [s^l and the clo- 
sure relation suggested by Rogers and Young [s^. The 
dimensionless Fourier transform of the pair correlation 
function c{k) = pc^ (fc) as a function of wave vector is 
plotted for different coupling constants F in Fig. [TJ 

Whereas the excess free energy in the RY approxima- 
tion to the DFT, Eq. ((TH]), requires, in general, the com- 
plete correlation function c{k), the respective function in 
the PFC models, Eq. ((^D)) . only needs the parameters 
Co, 6*2, and C4 as an input. The latter can in principle 
be obtained from c{k) in different ways. For our pur- 
pose they are chosen according to a series expansion of 
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c{k) in terms of k"^ about the correlation function's first 
maximum at fc = fc* up to second order in k^, i.e., 

c{k) ~Co + C2k^ + Cik'^ 

1 (32) 
= c(r) + {k^ - k*'^)c'{k*) + -(fc2 - k*yc"{k*) . 

Here, primes denote derivatives with respect to fc^. An- 
other way to determine the coefficients [25| is a a fit which 
reproduces the isothermal compressibility at fc — > 0, the 
bulk modulus, and the lattice constant of the crystal. 

As a subsequent motivation of the suggested fit and 
also for an estimate of the phase behavior, we calculated 
the Helmholtz free energy of the PFC2 model analytically 
in the one-mode approximation [SS^ . Within this approx- 
imation, the two-dimensional density field is assumed to 
be sinusoidal and hexagonally symmetric, i.e.. 



— cos ikx) 



cos 



\/3fca 



ky 



(33) 

with an amplitude A and a nearest neighbor distance of 
a = 27r/fc. Eq. ^ together with eqs. ([20]) and ([26]) yield 
the free energy per particle 



T{A,k) A^ 



N 



512 



{ 15^2 - IQA -)- 96 [1 - c(fc)] } . (34) 



Minimization with respect to k and A gives the equilib- 
rium wave number k = k* and the equilibrium amplitude 



A*[p) 



0, 



(20c(fc*;p) - 19) 



1/2 



1 



c{k*]p) < Cf 
c{k*;p) > Cf , 
(35) 

where c/ = 43/45 = 0.956. The first line corresponds 
to the stable fluid and the second to the stable crystal 
phase, respectively. For values of c„ < c{k*;p) < Cf, with 
c„ = 0.95, the crystalline density field is metastable, i.e., 
the free energy, Eq. (j34p , has a local, non-global minimum 
at a finite amplitude A. For values c{k*;p) < Cu, the 
crystal is unstable towards collapse. As the approximate 
free energy, Eq. in the PFC model, is governed by 

c(fc*), a proper representation of the latter is on order, 
and a series expansion of c{k) about k* appears natural. 

The RY approximation to the DFT predicts a stable 
crystal for F > F^^"^ « 36.5 and a metastable crystal 
for F„ < F < F/, with F„ « 31 [11]. The coupling 
constant at freezing F f corresponds to a maximum value 



DFT 

/ 



0.843, 



of the correlation function of c{k*;Tf) =: c 
which is substantially smaller than the value of c/ 
0.956, obtained in the one-mode approximation. On the 
other hand, extrapolation of c(fc*;F) to a value of c/ = 
0.956 yields a freezing transition within the one-mode 
approximation to the PFC2 model of Ff « 100, as can be 
seen in Fig.[2l In order to obtain a similar stability regime 
of the crystalline solution to the DDFT and to the PFC2 
and PFCl models, the excess free energy of the PFC 
models is rescaled by a factor of / = 1.15, independent 




FIG. 2: The maximum of the correlation function c{k*;T) 
versus the interaction strength F. The arrows at Fu and F/ 
bracket the range of metastability of the crystalline density 
field, obtained from the RY approximation to the DFT. 



of F. The phase behavior of both PFC models without 
the constraint on the functional form of the density field 
of Eq. is similar but slightly different than in the one- 
mode approximation; it is discussed in Subsection IVI Cl 



B. Non-equilibrium dynamics: Setup description 

In order to measure the propagation front velocities 
predicted by the DDFT, Eq. (fTT]) is numerically solved 
on a rectangular periodic box of a fine grid with ~ 64 
grid points per nearest-neighbor distance a. A finite dif- 
ference method with variable time step is applied. The 
convolution integrals are solved using the method of fast 
Fourier transform. The two versions of the PFC model, 
eqs. (|22l) and (|27|) . are solved by finite element methods 
with a variable timestep. The sixth order partial differ- 
ential equations are solved semi-implicitly as a system of 
second order partial differential equations 10']. 

We study the propagation dynamics of the front of a 
linear array along the y-direction, which is a cutout of 
a perfect hexagonal crystal and comprises, at t — 0, a 
number of s infinite rows of particles centered about x — 
0, as can be exemplarily seen from the density map for 
t = in Fig.[3]for F = 60. The number of crystalline rows 
s is chosen larger than the critical nucleus to guarantee 
crystal growth [2^ for F > F/ or reasonably large in 
order to study crystal shrinkage for F < F/|60']. The 
size of the periodic, rectangular box is therefore chosen 
integer multiples of the lattice spacing of the perfectly 
ordered hexagonal crystal. 



128(V3/2) a X a . 



(36) 
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The nearest-neighbor distance is fixed to its ideal value 



= (2/V3)i 



(37) 



which is very close to the equilibrium value of the one- 
mode approximation to the PFC2 model, a = 27r/fc*and 
to the equilibrium value of the RY approximation to the 
DFT for all values of T [H. The initial density field is 
given by 

p(r, i = 0) = [1 - hi\x\ - R)] pc(r) + hi\x\ - R)p . (38) 

Here, h{x) is a smoothed approximation to the Heaviside 
step function. Pc(r) is the infinite, stable or metastable, 
crystalline density field with constrained lattice constant 
a and with the [ll]-orientation parallel to the x-axis, 
which is symmetric about a: = 0. The size of the initial 
nucleus is defined by i? = (•\/3/4)sa, with s the num- 
ber of crystalline rows parallel to the y-axis. /i(|a;| — R) is 
therefore chosen to cut through valleys of the density field 
guaranteeing the overall particle number to be fixed, i.e., 
jy dr/9(r, t = {)) = p. For the number of crystalline 
rows for the different coupling constants, see |60| . 

The initial density field p(r, t — 0) can be thought of as 
an equilibrium density field with an appropriately chosen, 
though somewhat artificial, external potential V{y) p^ . 
An experimentally more feasible and thus more realistic 
setup has been suggested in reference |2^: The array of 
tagged particles was first, i.e., for times t < 0, held fixed 
in a thermodynamically stable, equilibrated fluid of den- 
sity p at a coupling constant F < well below freezing. 
For the equilibration of the fluid, Eq. (jl7p was numeri- 
cally solved fixing the tagged particles by deep parabolic 
external potentials at the tagged particle positions — in 
an experiment this could be achieved by using optical 
At time t = Q the external pinning poten- 



tweezers 

tial was turned off and, at the same time, the system was 
instantaneously quenched to a a higher coupling constant 
F. Experimentally, the instantaneous quench is easily 
achieved by increasing the homogeneous external mag- 
netic fleld B. Both protocols lead to different short-time 
dynamics [t < (pD)^^] due to the initial difference in the 
density field close to the incipient cluster. However, for 
longer times the density fields are indistinguishable (data 
not shown) thus leading to the same propagation front 
velocities which are of interest in this work. 



s 



1.5 



0.5 



■0.5 

-1 



pgio[p(r)/p] 



FIG. 3: (Color online) Snapshots of of the dimensionless, log- 
arithmic density field logj^Q[p(r, t)/p] of a linear nucleus of 
initially s = 5 (DDFT; top panel) or s = 11 (PFCl model; 
bottom panel) infinite rows of hexagonally crystalline parti- 
cles at coupling constant F = 60, as obtained from the DDFT 
(top panel) and the PFCl model (bottom panel). The upper 
and the lower four maps show the density fields each at times 
t/TB = 0, 0.5, 1, 1.5 (from top to bottom). For better visibility 
and exploiting the symmetry of the density field, the images 
display twice the right half of the system's central region of 
dimensions 35^/3/2a x 2a. 
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FIG. 4: The y-average density profile py(x,t — tb) obtained 
from the DDFT (top panel) and the PFCl model (bottom 
panel) for the same coupling constant F = 60 and the same 
initial conditions as in Fig. [3] at the time t/rs — 1. 



C. Non-equilibrium dynamics: Results 

In Fig.[3l we display the density field p{r, t) as obtained 
from the DDFT and from the rescaled PFC2 model for 
F = 60 at four different times t/rs = 0, 0.5, 1, 1.5, where 
we chose the Brownian time scale tb = {pD)~^ as the 
unit time scale. Since the density field is symmetric with 
respect to the a;-axis we concentrate on the region a; > 
in the following. Fig. [4] displays the corresponding y- 
averaged density profile, 

Py{x,t) = Ly^Uyp{v,t), (39) 



at a short time t — tb, which is large enough for Py{v, t) 
to be insensitive on the details of the same field at time 
t = {]. Density fields obtained within the rescaled PFCl 
model are qualitatively very similar to the ones of the 
PFC2 model; they are therefore not displayed in the 
present paper. 

It can be ascertained from Fig. [3] that after short time 
the propagating crystal fronts approximately establish 
steady states, which eventually change on larger time 
scales t ^ Tb '62*1 . As can also be seen from Figs. [3] 
and m the crystalline density fields behind the crystal 
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FIG. 5: (Color online) Propagation front velocities of a linear 
crystal front in the [ll]-direction, Vf{r), measured shortly af- 
ter the quench (see text), as a function of relative coupling 
constant AF = F — F/, with F/ the respective coupling con- 
stant of freezing, obtained within the DDFT (stars), the PFCl 
model (circles), and the PFC2 model (triangles). Lines are 
guides to the eye. Inset: The same velocities as a function of 
F. 



fronts, which are close to equilibrium in all three models, 
are modulated much stronger about the average density p 
in the DDFT than in the two rescaled PFC models. This 
difference goes along with an almost sinusoidal density 
field in the PFC models, approximately equal to the field 
assumed in the one-mode approximation, Eq. which 
is due to the truncation of the expansion of i^cx[/o] in fc^ 
[Eq. (HOI)]; this density field is in contrast to the one of 
the DDFT, which is approximately given by a superpo- 
sition of Gaussians centered about the lattice vectors of 
the hexagonal lattice. Another qualitative difference be- 
tween the two models concerns the crystal-melt interface: 
the width of the crystal front obtained within the DDFT 
is substantially smaller than in the PFC models; whereas 
the former is of the order of Aa; ^ 5p^^^^, the latter ap- 
proximately amounts Aa; ~ 25p~^/^. In both models, 
the widths are relatively insensitive towards changes in 
coupling constant F (data not shown). 

In order to quantify the crystal front propagation, the 
position of the diffuse front, Xf{t), is extracted as the 
maximum a;-position at which the density field exceeds 
the value 2p, in the DDFT, and as the inflection point 
of the envelope function to py{x,t) in the PFC model, 
respectively (data not shown). Within the DDFT, the 
front position is thus situated at a position, which is at 
the same time slightly larger than the one in the PFC 
models. However, this does not affect the crystal front ve- 
locity Vf{T) = dxf{t)/dt, which was measured as a func- 
tion of coupling constant F for the three different models 
under study at a short distance from the incipient front, 
Xf{t = 0) + 15/9-1/2 < xf{t) < Xf{t = 0) + 20p-i/2 (ggg 
Fig. O. The front velocity at long times is not consid- 
ered in this paper. Comparing the three different curves 



of Fig. [5l the following observations are made: 

(i) In all three models, w/(F) increases monotonically 
from with increasing difference in coupling constant 
AF = F — F/, where the freezing constant is Tf = 36.5 
within the DFT [H], F/ « 33 in the rescaled PFCl 
model, and F/ ft! 38 within the rescaled PFC2 model. 
The three coupling constants do not agree because the 
rescaling of the excess free energy was chosen to bring 
only the freezing constants of the more approximate one- 
mode approximation to the PFC2 model into agreement 
with the freezing constant of the DFT. The same scaling 
is then used for the PFCl and PFC2 models without the 
constraint of a sinusoidal density field, Eq. ([33]) . The dif- 
ference between the respective values within the PFCl 
and the PFC2 models is due to the difference in their 
respective ideal free energies -Fid[p], Eqs. (fl^ and pSj) . 
Moreover, as the ideal free energy functionals of the DFT 
and the PFCl model are equal, the freezing constant of 
the rescaled PFCl model is lower than the one of the 
rescaled PFC2 model, which is reminiscent of the smaller 
distance of the respective freezing constants within the 
non-rescaled models (obtained from an extrapolation of 
c(fc); data not shown). 

(ii) For negative AF, with F still within the metastabil- 
ity regime, F„ < F < F^, the propagation front velocity 
is negative and a retraction of an almost steady-state 
crystal front is observed. 

(iii) For AF > 0, the front velocities obtained within 
the two rescaled PFC models bracket the respective val- 
ues obtained within the DDFT, the PFCl model being 
sHghtly closer to the DDFT as the PFC2 model. All 
three models yield similar front velocities, although the 
corresponding steady-state density fields are quite differ- 
ent (see Fig. U). However, the non-monotonicity of the 
propagation front velocities w/(AF) for the same value 
of AF with increasing degree of approximation (from the 
DDFT via the PFCl model to the PFC2 model) points 
to a cancellation of errors in the approximate PFC mod- 
els. Still, as would have been expected, the velocities ob- 
tained from the less approximate PFCl model are closer 
to the results of the DDFT than those of the more ap- 
proximate PFC2 model. 



VII. CONCLUSIONS 

In conclusion, we have derived the phase field crys- 
tal (PFC) model from microscopic dynamical density 
functional theory (DDFT) appropriate for colloidal dis- 
persions, which are governed by completely overdamped 
Brownian dynamics. The ordinary phase field crystal 
model (called PFC2 model) arises from a constant mo- 
bility approximation and an expansion of the ideal gas 
entropy in terms of density. Both approximations can be 
avoided yielding a variant of the phase field crystal model 
(called PFCl model), which requires the same computa- 
tional effort as the PFC2 model. 

Comparing the two phase field crystal models to the 
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full solution of the dynamical density functional theory, 
agreement could only be obtained by an empirical scal- 
ing factor in the free energy. On the one hand, this im- 
plies that phase field crystal models have to be used with 
care and lack full ab initio precision. If a suitable scal- 
ing is accepted, there is good overall agreement in the 
growth velocity of a crystalline front, where the PFCl 
model yields slightly better agreement than the PFC2 
model. The corresponding density fields, however, differ 
vastly in terms of the sharpness of the crystalline peaks. 
We therefore conclude that the phase field crystal model 
gives a qualitative reliable description of the trends in 
crystal growth processes but lacks high precision. 

Future work should focus, first, on three-dimensional 
crystalline fronts and on the dynamics and annealing of 



crystalline defects, where both DDFT and PFC models 
should be compared as well. Second, the behavior of the 
setup at intermediate and long times deserves more de- 
tailed study. Finally, full Brownian dynamics computer 
simulations [24] are needed to obtain reference data for 
a test of the underlying adiabatic approximation in the 
DDFT. 
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